/*

Uses Vertex Coding to find proton pathways

Routine described by J. Scott, T. Ideker, R. M. Karp, and R. Sharan, 
"Efficient algorithms for detecting signalling pathways in protein interaction networks,"
Lecture notes in computer science, 3500, 1-13 (2005).
before colors are implemented

This is the same technique used in 
M. A. Gomez, D. Shepardson, L. T. Nguyen, and T. Kehinde, "Periodic long range proton conduction pathways in pseudo-cubic and orthorhombic perovskites",
Solid State Ionics, 213, 8 (2012)

Units:
Energies are in eV
Prefactors are in /ps.

Used in 
K.S. Gomez-Haibach and M.A. Gomez: Revised Centrality Measures Tell a Robust Story of Ion Conduction in Solids. 
The Journal of Physical Chemistry B 127, 9258 (2023).
and adapted by the Gomez Group at Mount Holyoke College for this tutorial 2025
*/

// Headers
#include<stdio.h>
#include<stdlib.h>
#include <math.h>

//Numerical Recipies

#include "nrutil.h"

//Errors
void MemoryError() {fprintf(stdout,"Error!! Ran out of Memory!\n");}

//routine to sort paths based on weight
void sortpaths(int pathsize, int **path,int numberpaths,double *weight,double *barrier,int *lengths,int *direction,double **distStore,int *BM)
{
    int i,j,k;
    int *temporarypath;
    temporarypath=ivector(1,pathsize);
    double temporaryweight,temporarybarrier,tempBM;
	int temporarylength,temporarydirection;
    double *temporaryDist;
    temporaryDist=dvector(1,3);

    for (j=2;j<=numberpaths;j++) {
		temporaryweight=weight[j];
		temporarybarrier=barrier[j]; 
		tempBM=BM[j];  
		temporarylength=lengths[j];
		temporarydirection=direction[j];
        for (k=1;k<=pathsize;k++) temporarypath[k]=path[j][k];
		for (k=1;k<=3;k++) temporaryDist[k]=distStore[j][k];
        i=j-1;
        while (i>0 && weight[i]>temporaryweight) {
			weight[i+1]=weight[i];
			barrier[i+1]=barrier[i];
			BM[i+1]=BM[i];
			lengths[i+1]=lengths[i];
			direction[i+1]=direction[i];
            for (k=1;k<=pathsize;k++) path[i+1][k]=path[i][k];
                for (k=1;k<=3;k++) distStore[i+1][k]=distStore[i][k];
                i--;
        }
        weight[i+1]=temporaryweight;
		barrier[i+1]=temporarybarrier;
		BM[i+1]=tempBM;
		lengths[i+1]=temporarylength;
		direction[i+1]=temporarydirection;
        for (k=1;k<=pathsize;k++) path[i+1][k]=temporarypath[k];
        for (k=1;k<=3;k++) distStore[i+1][k]=temporaryDist[k];
    }
     free_ivector(temporarypath,1,pathsize);
     free_dvector(temporaryDist,1,3);
}

int main()
{
	
//Reading in the number of binding sites 

int numberofsites;
int i,j,k,n,ii,mm;
fprintf(stdout,"\n********\nBEGINNING OF Vertex CODING	SETUP\n********\nHow many binding sites are there in this system?\n"); 
fscanf(stdin,"%d\n",&numberofsites);
fprintf(stdout,"%d\n",numberofsites); 

// Setting up Arrays  
//  Energy[i] will have the energy of the ith binding site 
//  TransitionState[i][j] will have the energy of the transition state between binding sites i and j
// RateConstant[i][j] will store the rate constants
//  Moves[i][j] will indicate the type of move between binding sites i and j
// It is equal to 0 for no move, 1 for rotations, 2 for transfers, 3 for inters

double *Energy;
Energy=dvector(1,numberofsites);
double **TransitionState;
TransitionState=dmatrix(1,numberofsites,1,numberofsites);
double **RateConstant;
RateConstant=dmatrix(1,numberofsites,1,numberofsites);
int **Moves;
Moves=imatrix(1,numberofsites,1,numberofsites);


// Setting up the files from which the Energy[i], Prefactor[i][j], TransitionState[i][j], and Moves[i][j] 
// arrays above will be read.  

FILE *TSFile;
TSFile =fopen("Energies.TS","r"); 
if (TSFile ==0)
  {
    fprintf(stderr,"Can't open File Energies.TS");
    exit(1);
  }

FILE *EnergyFile;
EnergyFile=fopen("Energies.binding","r");
if (EnergyFile==0)
  {
          fprintf(stderr,"Can't open File Energies.binding");
          exit(1);
   }

//Setting up the temperature and beta
//1.3806503e-23 is Boltzmann's constant in units of J/K
//6.24150974e18 is the conversion from J to eV
double beta, temperature, kb;
kb=1.3806503e-23*6.24150974e18;
fprintf(stdout,"What temperature would you like to run at? (K)\n");
fscanf(stdin,"%lf\n",&temperature);
fprintf(stdout,"%lf Kelvin\n",temperature);
beta=1/(kb*temperature);
fprintf(stdout,"beta is %lf in units of 1/eV\n",beta);

//Reading binding and transition state energies, prefactors, and move type files
//Calculating the rate constant matrix

double **Bpositions;
Bpositions=dmatrix(1,numberofsites,1,3);

int isite;


fprintf(stdout,"Reading binding site positions and energies\n");
for(i=1;i<=numberofsites;i++)
{
 fscanf(EnergyFile,"%d %lf %lf %lf %lf\n",&isite,&Bpositions[i][1],&Bpositions[i][2],&Bpositions[i][3],&Energy[i]);
 //fprintf(stdout,"Site %d %lf %lf %lf Energy %lf\n",isite,Bpositions[i][1],Bpositions[i][2],Bpositions[i][3],Energy[i]);
 for(j=1;j<=numberofsites;j++)
   {
    TransitionState[i][j]=0.0;
    Moves[i][j]=0;
    RateConstant[i][j]=0.0;
   }
}
//now the ij files
int numberTS,numberRot,numberIntra,numberInter,*from,*to;
double TSenergy;
int movesEntry;
//fscanf(TSFile,"%d\n",&numberTS);
numberTS=216;
from=ivector(1,numberTS); to=ivector(1,numberTS);
numberRot=0; numberIntra=0; numberInter=0;
fprintf(stdout,"Reading TS site, move type, energies\n");
for (i=1;i<=numberTS;i++) 
{
 fscanf(TSFile,"%d %d %d %lf \n",&from[i],&to[i],&movesEntry,&TSenergy);//reads transition state energies
 Moves[from[i]][to[i]]=movesEntry;
 TransitionState[from[i]][to[i]]=TSenergy;
 //fprintf(stdout,"%d %d %d %lf\n",from[i],to[i],Moves[from[i]][to[i]],TransitionState[from[i]][to[i]]);
 TransitionState[to[i]][from[i]]=TransitionState[from[i]][to[i]]; 
 Moves[to[i]][from[i]]=Moves[from[i]][to[i]]; 

 if (Moves[to[i]][from[i]]==1) {
		 numberRot=numberRot+1;
 } else {
	 if (Moves[to[i]][from[i]]==2) {
		numberIntra=numberIntra+1;
	 } else {
		if (Moves[to[i]][from[i]]==3) {
				numberInter=numberInter+1;
		} else {
			fprintf(stdout,"There is an unaccounted move\n");
			exit(1);
		}
	 }
 }
}
 
 fprintf(stdout,"There are %d intra, %d rot, %d inter moves\n",numberIntra,numberRot,numberInter);

//Finding the lowest energy site
double minEnergy;  minEnergy=0.0;//large number as all minima have negative energies
int imin;
for (i=1;i<=numberofsites;i++) {
        if (minEnergy>Energy[i]) {
			minEnergy=Energy[i];
			imin=i;
		}
}
fprintf(stdout,"Mininum Energy is %lf and occurs for binding site %d\n",minEnergy,imin);

double boltznorm; boltznorm=0.0;
//Boltmann normalization factor
for (i=1;i<=numberofsites;i++) boltznorm=boltznorm+exp(-beta*(Energy[i]-minEnergy)); 
// if we had frequencies, each term would need to be divided by the product of frequencies at the ith site
//fprintf(stdout,"boltz norm %lf\n",boltznorm);

//fprintf(stdout,"Site Position Energy Probability\n");
//for (i=1;i<=numberofsites;i++) 
	//fprintf(stdout,"%d %lf %lf %lf %lf %e\n",i,Bpositions[i][1],Bpositions[i][2],Bpositions[i][3],Energy[i]-minEnergy,exp(-beta*(Energy[i]-minEnergy))/boltznorm);
 // probability would also need to be divided by the product of frequencies at each state if we had frequencies

//Calculating rate constant matrix
fprintf(stdout,"Rate constants and prefactors\n");
//fprintf(stdout,"from to MoveType TS energy  for/back barrier forward backward ktst  forward backward A\n");

for (i=1;i<=numberofsites;i++) {
        for (j=i+1;j<=numberofsites;j++) {
                if (Moves[i][j]>0) {
                        RateConstant[i][j]=exp(-beta*(TransitionState[i][j]-Energy[i]));//would need to multiplied by a prefactor (ratio of frequencies at minimum to those at TS if we had frequencies
                        RateConstant[j][i]=exp(-beta*(TransitionState[j][i]-Energy[j]));
                        //fprintf(stdout,"%d %d %d %lf %lf %lf %lf %lf\n",i,j,Moves[i][j],TransitionState[j][i]-minEnergy,TransitionState[i][j]-Energy[i],TransitionState[j][i]-Energy[j],RateConstant[i][j],RateConstant[j][i]);
                }
        }
}

//Calculating weights from rate constants
// weight=-ln(probabilities)

double **weight;
weight=dmatrix(1,numberofsites,1,numberofsites);

//probabilities
double sum,minimum;
minimum=100000.0;  //very large number
double **probabilities;
probabilities=dmatrix(1,numberofsites,1,numberofsites);

//fprintf(stdout,"\ni j Probabilities[i][j] RateConstant[i][j]\n");
for(i=1;i<=numberofsites;i++){
    sum=0.0;
    for (j=1;j<=numberofsites;j++) sum=sum+RateConstant[i][j];
    for(j=1;j<=numberofsites;j++){
    	probabilities[i][j]=RateConstant[i][j]/sum;
		if (probabilities[i][j]>0) {
			weight[i][j]=-log(probabilities[i][j]);
		} else {
			// 0 is a flag for infinite weight
			weight[i][j]=0;
		}
     	//if (probabilities[i][j]>0) fprintf(stdout,"%d %d %e %e\n",i,j,probabilities[i][j],RateConstant[i][j]);
    }
}

// input and start
int initialBindingSite,fromSite,toSite,lengthofpath,length1,length2;
fprintf(stdout,"What is first binding site number to consider?\n");
fscanf(stdin,"%d\n",&fromSite);
fprintf(stdout,"%d\n",fromSite);
fprintf(stdout,"What is the final binding site number to consider?\n");
fscanf(stdin,"%d\n",&toSite);
fprintf(stdout,"%d\n",toSite);
fprintf(stdout,"What is the first desired path length?\n");
fscanf(stdin,"%d\n",&length1);
fprintf(stdout,"%d\n",length1);
fprintf(stdout,"What is the final desired path length?\n");
fscanf(stdin,"%d\n",&length2);
fprintf(stdout,"%d\n",length2);
double criticalDistance;
fprintf(stdout,"What is the critical distance?\n");
fscanf(stdin,"%lf\n",&criticalDistance);
fprintf(stdout,"%lf\n",criticalDistance);

int maxNumberConnections;
maxNumberConnections=5; //2 rotations, 2 intra, and 1 or 0 inter 
int MaxNumberPossibilities;
MaxNumberPossibilities=numberofsites;
for (i=2;i<=length2;i++) MaxNumberPossibilities=MaxNumberPossibilities*maxNumberConnections;
fprintf(stdout,"MaxNumberPossibilities %d\n",MaxNumberPossibilities);
//reset since most of the time we don't realize as many as calculated
MaxNumberPossibilities=10000000;
fprintf(stdout,"Reset MaxNumberPossibilities %d\n",MaxNumberPossibilities);
	
int **possibilities;
possibilities=imatrix(1,MaxNumberPossibilities,1,length2);
double *temporaryWeight;
temporaryWeight=dvector(1,MaxNumberPossibilities);

int numberPossibilities,nbestpath;
int *siteSet;	
siteSet=ivector(1,numberofsites);

bool deadEnd;
int flagSite;
int *newPossibilities;
newPossibilities=ivector(1,maxNumberConnections);

int countNewPossibleSites;

double *totalweight;
totalweight=dvector(1,MaxNumberPossibilities);
double *maxBarrier;
maxBarrier=dvector(1,MaxNumberPossibilities);
int *direction;
direction=ivector(1,MaxNumberPossibilities);
int **paths;
paths=imatrix(1,MaxNumberPossibilities,1,length2);
int numberofpathsfound;
numberofpathsfound=0;
int *lengthPaths;
lengthPaths=ivector(1,MaxNumberPossibilities);
for (i=1;i<=MaxNumberPossibilities;i++) lengthPaths[i]=0;
double average,normalizing,property,temp;
average=0.0; normalizing=0.0;
double avgX,avgY,avgZ,normX,normY,normZ;
avgX=0.0; avgY=0.0; avgZ=0.0;
normX=0.0; normY=0.0; normZ=0.0;

double distTemp,rdist[4]; //double pathDistance;
double **distStore;
distStore=dmatrix(1,MaxNumberPossibilities,1,3);
int flagMove;
double nT,nI,nR;
nT=0.0; nI=0.0; nR=0.0;
int *BM;
BM=ivector(1,MaxNumberPossibilities);

//Vertex Coding Starts
sum=0.0;
//start loop for all sites and lengths
  for (initialBindingSite=fromSite;initialBindingSite<=toSite;initialBindingSite++) {
    for (lengthofpath=length1;lengthofpath<=length2;lengthofpath++) {
	//initiallize number of possiblities for this specific initial site and length of path
	numberPossibilities=1;
	//initialize all paths to start at the initialBindingSite chosen
	// and with zero weight
	for (i=1;i<=MaxNumberPossibilities;i++) { 
		possibilities[i][1]=initialBindingSite;
  		temporaryWeight[i]=0.0;
	        for (j=2;j<=lengthofpath;j++) possibilities[i][j]=0;
	}
	//loop throught all possiblities - numberPossiblities starts at 1 but gets advanced in loop as paths are made
	for (n=1;n<=numberPossibilities;n++) {
	   //initializing set of possible sites
	    for (i=1;i<=numberofsites;i++) siteSet[i]=i; 	
	   //initializing ii and advancing it to first empty spot in nth possibility while removing used sites from site set
	    ii=1;
	    while (possibilities[n][ii]>0) {
			siteSet[possibilities[n][ii]]=0;//removing site 
			ii=ii+1; //advancing ii 
	    }
	    //Building of paths or looping through the stages
	    for (i=ii;i<=lengthofpath;i++) {  
		countNewPossibleSites=0; //initializing count of new possible sites 
		deadEnd=true; //initializing flag for dead end mark - changed only if continuation paths are found 
		//Looping through sites to see all possible next steps
		for (j=1;j<=numberofsites;j++) {
			if (weight[possibilities[n][i-1]][j]>0) {  //checking if new possible site is connected 
				flagSite=0; //initializing flag on site availability - will change if site is available 
				for (k=1;k<=numberofsites;k++) {
				    if (j==siteSet[k]) {flagSite=j; k=numberofsites+1;}  //site is available, exit loop 
			        }
				if (flagSite>0) {  //flag indicates that site was available 
					countNewPossibleSites=countNewPossibleSites+1; //increase count of possible sites
					newPossibilities[countNewPossibleSites]=flagSite; //keep site in newPoss vector 
					deadEnd=false;  //signals that this is not a deadend i.e a site was available
				} 
			}
		}
		//now put newpossibilities into possibilities matrix
		//as long as it is possible to continue path 
		if (deadEnd==false) {
		    //new possibility is accepted
		    possibilities[n][i]=newPossibilities[1];
		   //remove vertex from siteSet
		    siteSet[possibilities[n][i]]=0;  //zero is the marker for unavailable site 
		    for (k=2;k<=countNewPossibleSites;k++) {  // putting additional possibilities in matrix 
			for (mm=1;mm<=i-1;mm++) possibilities[numberPossibilities+k-1][mm]=possibilities[n][mm];//old sites
			possibilities[numberPossibilities+k-1][i]=newPossibilities[k]; //new site 
		    	temporaryWeight[numberPossibilities+k-1]=temporaryWeight[n]+
				weight[possibilities[numberPossibilities+k-1][i-1]][possibilities[numberPossibilities+k-1][i]];
		    } 
		    numberPossibilities=numberPossibilities+countNewPossibleSites-1; //adjustment of upper n bound 
		    temporaryWeight[n]=temporaryWeight[n]+weight[possibilities[n][i-1]][possibilities[n][i]];
		    //fprintf(stdout,"Updating number of possibilities to %d\n",numberPossibilities);
		 // }

		 } else {
			//deadEnd was false - exit path
			i=lengthofpath+1;
		 }
	      }//end of i loop for making path 
        } // end of n loop that went over all possiblities 

//Just storing paths that are span system now
//and start and end at the same point
          
   for (n=1;n<=numberPossibilities;n++) {
	if (possibilities[n][lengthofpath]>0) { //check that it's long enough

	   if (weight[possibilities[n][lengthofpath]][initialBindingSite]>0) {  //check that it connects to original 

 	  		//pathDistance=0.0;
          		for (k=1;k<=3;k++) rdist[k]=0.0;
          		//adding vector x, y, and z displacements over path
          		//taking into account periodicic boundary conditions on individual segments
          		// this way only for loops do all three components go to zero
          		//for things spanning cell, one component goes to 1.0
          		for (i=1;i<=lengthofpath-1;i++) {
              			for (k=1;k<=3;k++) {
                			distTemp=Bpositions[possibilities[n][i+1]][k]-
                        			Bpositions[possibilities[n][i]][k];
                			//periodicities in each section of path
                			if (distTemp<-0.5) distTemp=distTemp+1.0;
                			if (distTemp>0.5) distTemp=distTemp-1.0;
                			rdist[k]=rdist[k]+distTemp;
              			}
          		}
              		for (k=1;k<=3;k++) {
                		distTemp=Bpositions[possibilities[n][1]][k]-
                        		Bpositions[possibilities[n][lengthofpath]][k];
                		//periodicities in each section of path
                		if (distTemp<-0.5) distTemp=distTemp+1.0;
                		if (distTemp>0.5) distTemp=distTemp-1.0;
                		rdist[k]=rdist[k]+distTemp;
              		}
          		//pathDistance=sqrt(rdist[1]*rdist[1]+rdist[2]*rdist[2]+rdist[3]*rdist[3]);

          	 //if (pathDistance>criticalDistance) { // check that it's long enough 
          	 if ((abs(rdist[1])>criticalDistance)||(abs(rdist[2])>criticalDistance)||(abs(rdist[3])>criticalDistance)) { // check that it's long enough 

		   	property=0.0;   
		   	for (i=2;i<=lengthofpath;i++) {
			  temp=TransitionState[possibilities[n][i-1]][possibilities[n][i]]-Energy[possibilities[n][i-1]]; 
			  if (temp>property) {
				property=temp; 
			    flagMove=Moves[possibilities[n][i-1]][possibilities[n][i]]; 	
			  }
		    }
			temp=TransitionState[possibilities[n][lengthofpath]][initialBindingSite]-Energy[possibilities[n][lengthofpath]]; 
			if (temp>property) {
				property=temp; 
			        flagMove=Moves[possibilities[n][lengthofpath]][initialBindingSite]; 	
			}

			temporaryWeight[n]=temporaryWeight[n]+weight[possibilities[n][lengthofpath]][initialBindingSite];
			//ones to consider
		   	numberofpathsfound=numberofpathsfound+1;
			lengthPaths[numberofpathsfound]=lengthofpath;
		   	for (i=1;i<=lengthofpath;i++) paths[numberofpathsfound][i]=possibilities[n][i];
		   	totalweight[numberofpathsfound]=beta*(Energy[initialBindingSite]-minEnergy)+temporaryWeight[n]+log(boltznorm);
			//the +log(boltznorm) factor cancels out in the overall normalization
			//but we'll leave it here
		   	if (totalweight[numberofpathsfound]<minimum) {
				minimum=totalweight[numberofpathsfound]; 
				nbestpath=numberofpathsfound;
			}
			maxBarrier[numberofpathsfound]=property; 
			BM[numberofpathsfound]=flagMove;
			if (flagMove==1) nR=nR+exp(-totalweight[numberofpathsfound]);
			if (flagMove==2) nT=nT+exp(-totalweight[numberofpathsfound]);
			if (flagMove==3) nI=nI+exp(-totalweight[numberofpathsfound]);
			// averaging over all paths
			// that are long enought and connect back 
			average=average+maxBarrier[numberofpathsfound]*exp(-totalweight[numberofpathsfound]); 
			normalizing=normalizing+exp(-totalweight[numberofpathsfound]);
                        distStore[numberofpathsfound][1]=rdist[1];
                        distStore[numberofpathsfound][2]=rdist[2];
                        distStore[numberofpathsfound][3]=rdist[3];
   			if (abs(rdist[1])>criticalDistance) {
                          avgX=avgX+maxBarrier[numberofpathsfound]*exp(-totalweight[numberofpathsfound]);
                          normX=normX+exp(-totalweight[numberofpathsfound]);
			  direction[numberofpathsfound]=1;
                	}
                	if (abs(rdist[2])>criticalDistance) {
                          avgY=avgY+maxBarrier[numberofpathsfound]*exp(-totalweight[numberofpathsfound]);
                          normY=normY+exp(-totalweight[numberofpathsfound]);
						  direction[numberofpathsfound]=2;
                	}
                	if (abs(rdist[3])>criticalDistance) {
                          avgZ=avgZ+maxBarrier[numberofpathsfound]*exp(-totalweight[numberofpathsfound]);
                          normZ=normZ+exp(-totalweight[numberofpathsfound]);
						  direction[numberofpathsfound]=3;
                	}
			//fprintf(stdout,"numberofpathsfound=%d length=%d weight=%lf %d %d\n",numberofpathsfound,lengthPaths[numberofpathsfound],totalweight[numberofpathsfound],possibilities[n][lengthofpath],initialBindingSite); 
	        }  //end of pathDistance check 

   	   } //end of connecting back to original check 
	} //end of checking that the path did not have a dead end i.e. is less than expected path length 
    } //end of n loop over possiblities found 
 //fprintf(stdout,"%d numberofpathsfound=%d\n",initialBindingSite,numberofpathsfound);
  } // end of loop over path lengths 
} // end of loop over initial binding sites 
//Results
	fprintf(stdout,"\naverage limiting barrier=%lf\n",average/normalizing);
	fprintf(stdout,"fractionInXDirection=%lf fractionInYDirection=%lf fractionInZDirection=%lf\n",normX/normalizing,normY/normalizing,normZ/normalizing);
	fprintf(stdout,"\nMax barrier X is  %lf\nMaxBarrier Y is %lf\nMax Barrier Z is %lf\n",avgX/normX,avgY/normY,avgZ/normZ);
	fprintf(stdout,"nR=%lf nT=%lf nI=%lf Sum=%lf\n",nR/normalizing,nT/normalizing,nI/normalizing,(nR+nT+nI)/normalizing);
	fprintf(stdout,"There are %d such paths\n",numberofpathsfound);

       fprintf(stdout,"Best one has weight %lf \n",minimum);
       fprintf(stdout,"Length of best path is %d\n",lengthPaths[nbestpath]);
       n=nbestpath;
           for (i=1;i<=lengthPaths[n];i++) fprintf(stdout,"%d ",paths[n][i]);
           fprintf(stdout,"%d %lf ",paths[n][1],totalweight[n]);
           for (i=2;i<=lengthPaths[n];i++) {
                if (Moves[paths[n][i-1]][paths[n][i]]==1) fprintf(stdout,"R");
                if (Moves[paths[n][i-1]][paths[n][i]]==2) fprintf(stdout,"T");
                if (Moves[paths[n][i-1]][paths[n][i]]==3) fprintf(stdout,"I");
                if (Moves[paths[n][i-1]][paths[n][i]]==0) { fprintf(stdout,"error"); exit(1);}
           }
          if (Moves[paths[n][lengthPaths[n]]][paths[n][1]]==1) fprintf(stdout,"R");
          if (Moves[paths[n][lengthPaths[n]]][paths[n][1]]==2) fprintf(stdout,"T");
          if (Moves[paths[n][lengthPaths[n]]][paths[n][1]]==3) fprintf(stdout,"I");
          if (Moves[paths[n][lengthPaths[n]]][paths[n][1]]==0) {fprintf(stdout,"error"); exit(1);}


//for long paths - dont' sort just top here
       fprintf(stdout,"\nTo see sorted paths, comment out exit in program below this print line in code\n");
       exit(1); 
       fprintf(stdout,"\nNow the sorted paths are:\n");

//sorting through rest all paths
	sortpaths(length2,paths,numberofpathsfound,totalweight,maxBarrier,lengthPaths,direction,distStore,BM);

//printing out sorted paths
	double expectedTime;
	double eTavg; eTavg=0.0;
	sum=0.0;	
	for (n=1;n<=numberofpathsfound;n++) {
           sum=sum+exp(-totalweight[n])/normalizing;
	   if (exp(-totalweight[n])/exp(-totalweight[1])<0.1) fprintf(stdout,"Insignificant ");
		//insignificant means that the probabilty is 1/10 of the best
	   expectedTime=0.0;
	   for (i=1;i<=lengthPaths[n]-1;i++) 
		expectedTime=expectedTime+1/RateConstant[paths[n][i]][paths[n][i+1]];  
	   expectedTime=expectedTime+1/RateConstant[paths[n][lengthPaths[n]]][paths[n][1]]; 
	   eTavg=eTavg+expectedTime*exp(-totalweight[n]);
	   fprintf(stdout,"n=%d length=%d Path ",n,lengthPaths[n]);
           for (i=1;i<=lengthPaths[n];i++) 
		fprintf(stdout,"%d ",paths[n][i]); 
//will put in a weight that does not consider initial point for comparison
           fprintf(stdout,"%d weight=%lf(%lf) prob=%lf probsum=%lf ET=%lf ",paths[n][1],totalweight[n],totalweight[n]-beta*(Energy[paths[n][1]]-minEnergy)-log(boltznorm),exp(-totalweight[n])/normalizing,sum,expectedTime);

	   for (i=2;i<=lengthPaths[n];i++) {
		if (Moves[paths[n][i-1]][paths[n][i]]==1) fprintf(stdout,"R"); 	
		if (Moves[paths[n][i-1]][paths[n][i]]==2) fprintf(stdout,"T"); 	
		if (Moves[paths[n][i-1]][paths[n][i]]==3) fprintf(stdout,"I"); 	
		if (Moves[paths[n][i-1]][paths[n][i]]==0) { fprintf(stdout,"error"); exit(1);}	
           }
	  if (Moves[paths[n][lengthPaths[n]]][paths[n][1]]==1) fprintf(stdout,"R"); 	
	  if (Moves[paths[n][lengthPaths[n]]][paths[n][1]]==2) fprintf(stdout,"T"); 	
	  if (Moves[paths[n][lengthPaths[n]]][paths[n][1]]==3) fprintf(stdout,"I"); 	
	  if (Moves[paths[n][lengthPaths[n]]][paths[n][1]]==0) {fprintf(stdout,"error"); exit(1);} 	

	  if (BM[n]==1) fprintf(stdout," LimitingStep R"); 
	  if (BM[n]==2) fprintf(stdout," LimitingStep T"); 
	  if (BM[n]==3) fprintf(stdout," LimitingStep I"); 
      fprintf(stdout," Barrier1=%lf",maxBarrier[n]); 

	  fprintf(stdout," Disp=(%lf,%lf,%lf)",distStore[n][1],distStore[n][2],distStore[n][3]);

	  if (direction[n]==1) fprintf(stdout,"x-direction\n");
	  if (direction[n]==2) fprintf(stdout,"y-direction\n");
	  if (direction[n]==3) fprintf(stdout,"z-direction\n");
                fprintf(stdout,"H Hnumber x y z FE_TS BE_TS ETforstep\n");
                for (i=1;i<=lengthPaths[n]-1;i++) fprintf(stdout,"H H%d %lf %lf %lf %lf %lf %lf\n",
                        paths[n][i],
                        Bpositions[paths[n][i]][1],Bpositions[paths[n][i]][2],
                        Bpositions[paths[n][i]][3],
                        TransitionState[paths[n][i]][paths[n][i+1]]-Energy[paths[n][i]],
                        TransitionState[paths[n][i]][paths[n][i+1]]-Energy[paths[n][i+1]], 
						1/RateConstant[paths[n][i]][paths[n][i+1]]);
						i=lengthPaths[n]; 
				fprintf(stdout,"H H%d %lf %lf %lf %lf %lf %lf\n",
                        paths[n][i],
                        Bpositions[paths[n][i]][1],Bpositions[paths[n][i]][2],
                        Bpositions[paths[n][i]][3],
                        TransitionState[paths[n][i]][paths[n][1]]-Energy[paths[n][i]],
                        TransitionState[paths[n][i]][paths[n][1]]-Energy[paths[n][1]],
						1/RateConstant[paths[n][i]][paths[n][1]]);
						fprintf(stdout,"H H%d %lf %lf %lf\n",paths[n][1],
                        Bpositions[paths[n][1]][1],Bpositions[paths[n][1]][2],
                        Bpositions[paths[n][1]][3]);
	}
	//above printed positions of points in path for quick imaging but for paper images overal vesta files that include backbone optimization
//re-iterating average results
	fprintf(stdout,"\naverage=%lf \n",average/normalizing);
	fprintf(stdout,"frac_x=%lf frac_y=%lf frac_z=%lf\n",normX/normalizing,normY/normalizing,normZ/normalizing);
	fprintf(stdout,"\nMax barrier X is  %lf\nMaxBarrier Y is %lf\nMax Barrier Z is %lf\n",avgX/normX,avgY/normY,avgZ/normZ);
	fprintf(stdout,"nR=%lf nT=%lf nI=%lf Sum=%lf\n",nR/normalizing,nT/normalizing,nI/normalizing,(nR+nT+nI)/normalizing);
	fprintf(stdout,"ExpectedTimeAvg=%lf\n",eTavg/normalizing);
}
